function fitting()
% fitting.m - The purpose of this is to take an excel spreadsheet formatted
% the correct
%     way and from the data in this calculate TEq, DeltaHEq, DeltaGInact and DeltaGCat
%     to get a best fit (using least squares method) to the equation in calc.m. This is
%     a non-linear fit so uses a iterative process to get closer and closer to the best
%     fit and exit when it reaches near this.
%
%     James McDowall - jcmcdowall@gmail.com
%
%  Changed in July 2006 to output standard deviations  and confidence limits by Karin Bryan and
%  Nigel Goodhue, U of Waikato.
%
%  Changed in Sept 2006 to remove most of the for loops to speed it up by
%  Karin Bryan
%
%  Changed in Dec 2006 to incorperate t = 0 calculations and addition of a
%  meshgrid for the rates by Daniel Hailstone, UoW.

% This sets the constants
R = 8.31447215;
Kb = 1.38065812e-23;
h = 6.626075540e-34;

changeCutoff = .000001;
webpage = 1;      %This is a variable to determine if the program is embedded
                  %in a webpage, if it is there is no graphing and a better
                  %formatted output. 1 = not webpage, this functionality
                  %has not been totally developed at the moment but is here
                  %as a skeleton idea of what can be altered if this was to
                  %be used in a webpage as was a propsed idea..

clear file;
clear ext;
clear ZAct;

% Gets the user to select the file that they want to access, and edit.
[file ext] = uigetfile('.xls');
if isequal(ext,0) || isequal(file,0)
    return;
end
filename = [ext file];

% Starts the timer
tic;

%Gets the data from the file that the user inputted as the filename
%This file provides the data, plus initial guesses of the constant, and 
%and error tolerance.
[totalData allData E0 DeltaGCat DeltaGInact DeltaHEq TEq err highres]  = getData(filename);


%Determine if the system will preform a t = 0 calculation or not.
if (size(totalData,1) == 2) & (totalData(2,1) == 0)
    'Performing a time = 0 calculation.'
    t0calc = 1;
else
    t0calc = 0;
end


%Sets the Least Mean Square as the minimum wanted error, so that the loop
%will be entered into
LMS = err*1;

% This sets up the initial values for some of the variables
up = 0;                         %Allows the program to know which direction it is moving


addition = [0 0 0 0];
additionSize = size(addition);  %Size of the addition array (for loop counts)
additionSize = additionSize(2); %Changes to be equal to just the number of columns.
totDataSize = size(totalData);  %Size of Data (for loop counts)
lines = 4;                      %The number of lines to access from totalData
prevLines = 0;                  %Variable to store prevoius number of lines included (used in calculating percentage)

% This is a loop which gets a fit to the data, on each loop it includes
% more of the data in the fitting procedure each time, the reason for this
% is that while using a small amount of data the fitting happens quickly,
% so by using less data you can quickly get near the required values, then
% introduce more data and end up with an accurate fit relatively quickly
while true
    
    firstChange = .1;           %This variable is used so the loop can know when enough reductions in change have been made
    change = firstChange;       %Change is the variable that lets the loop know how much to step by each time

    clear Data;

    %Maps the temperature data to a new array
        Data(1,1:totDataSize(2)) = totalData(1,1:totDataSize(2));
        
      
    %This calculates the number of lines of data that will need to be jumped
    %betweeen the ones included in Data, to allow 'lines' lines to be
    %included
    jump = floor(totDataSize(1)/lines) + 1;


    %This is the algorithm that steps through the input data and pulls out
    %the lines that area jump apart, meaning that in Data at the end there
    %is onle 'lines' rows.
        z = 2:jump:totDataSize;
        y = 1:totDataSize(2);
        Data((z-2)/jump + 2,y) = totalData(z,y);
        

    dataSize = size(Data);          %Since the size of data changes each loop need to be recalculated
    
    value = [DeltaGCat DeltaGInact DeltaHEq TEq LMS*2]; %Restarts the value array, which
                                                        %contains the variable values and
                                                     %relavent mean square data in the
                                                        %5th column
    count = 0;                          %A variable to count the number of elements included in the loop
    value(1,5) = 0;             %The last column contains the mean square data


    %This finds the root mean square error for the top row of the value
    %array (same algorithm as described later).
    
    Temp = Data(1,2:dataSize(2));          %All temperature values
    Time = Data(2:dataSize(1),1);           %All time values        
    [TT,tt]=meshgrid(Temp,Time);            %Make a grid of all possibilities
    T=reshape(TT,(dataSize(2)-1)*(dataSize(1)-1),1);         %reshape grid into vector
    t=reshape(tt,(dataSize(2)-1)*(dataSize(1)-1),1);
    y=reshape(Data(2:dataSize(1),2:dataSize(2)),(dataSize(2)-1)*(dataSize(1)-1),1);
    
                                                
                                               
       %Finds the square of the difference between the datapoint and the equation with current variables   
      if t0calc == 1
             value(1,5) =sqrt(mean((calcT0(T,t,value(1,1),value(1,2),value(1,3),value(1,4),E0,h,Kb,R)-y).^2));   
      else
            value(1,5) =sqrt(mean((calc(T,t,value(1,1),value(1,2),value(1,3),value(1,4),E0,h,Kb,R)-y).^2));
      end
  
       
                    
     
    error = value(1,5);                 %The error is the mean square differences

    %This is the loop that finds the best fit to the data that is currently
    %in the array 'Data'. As it loops through the step size and position
    %changes, this loops is exited only when the step size has got to a
    %small enough size, or the error changes by a very small amount each time.

    while true
        
        %Sets addition array to inital values
        addition = [ 0 0 0 -1];
                  
        j = 2;                          %Variable to keep track of which line of value array is
                                        %currently being worked on

        %This loops steps through the possible next positions if only one
        %variable at a time can be altered, and each variable can be
        %altered only by the change size.
        while true;

            %For the current addition setting, this sets what each
            %variable would be
                value(j,1:4) = value(1,1:4) + addition.* change .* value(1,1:4);

            %This increments the addition array appropriately, only one
            %element at a time can not be 0, and they can only be negative
            %or postive 1. Changes through all the permutations
            %sequentially, once complete sets the 4th element to 3
            for i = additionSize:-1:1
                if addition(i) == -1
                    addition(i) = addition(i) + 2;
                elseif addition(i) == 1
                    if i == 1
                        addition(additionSize) = 3;
                    else
                        addition(i) = 0;                      
                        addition(i-1) = -1;
                        break;
                    end
                end
            end

            value(j,5) = 0;                 %Initialse to 0, so it can be added to.

            %This calculates the square difference between the all the
            %datapoints and the equation with the current coefficients
            %and adds them all up, they are then divided by the number of
            %datapoints and the square root of that is taken, to give the
            %root mean square error.
           if  t0calc == 1
            value(j,5) =sqrt(mean((calcT0(T,t,value(j,1),value(j,2),value(j,3),value(j,4),E0,h,Kb,R)-y).^2)); 
           else
            value(j,5) =sqrt(mean((calc(T,t,value(j,1),value(j,2),value(j,3),value(j,4),E0,h,Kb,R)-y).^2));   
           end
                                                %Finds the square of the difference
                                                %between the datapoint and the
                                                %equation with current variables        
           % error = value(j,5);               %The error is the mean square differences
            value(j,5) = value(j,5)/length(y);

            %This checks if all the permutations of addition have been used
            if addition(additionSize)==3
                break;
            end

            j = j + 1;                      %Moves to next line of the array

  
        end
        [number pos] = min(value(:,5)) ;    %Finds the coefficient values with
                                            %the smallest mean square
                                            %difference, which is best fit

        value(1,:) = value(pos,:);          %Moves the best fit to the top
        %                                    %where it will be the start for
        %                                    %the next iteration
        error2 = error;
        error = value(1,5);
  

        %This checks if the best fit is actually the top row, if it is then
        %the size of change is changed to see if there is a better fit in a
        %different region. First moves to a larger change because when
        %stepping a larger step is better to reduce number of iterations,
        %if there is no better fit in a larger step size then step is
        %reduced to below initial value.
  
        if pos == 1 || value(pos,5) == value(1,5) %if the new value are as good as or equal to the last ones
            %This is if the previous move was an increase
            if up == 1
                change = change * 4/9;      %Reduces the step size
                up = 0;                     %Resets the step, it last went down
            %This is if the previous move was not an increase
            else
                change = change * 1.5;      %Increases the step size
                up = 1;                     %The last step was now a increase
            end
            %This checks to see if the change in error is too small, if it is
            %then the step size is increased.
        elseif (abs(error - error2))/error < changeCutoff
            %This is if the previous move was an increase
            if up == 1
                change = change * 4/9;      %Reduces the step size
                up = 0;                     %Resets the step, it last went down
            %This is if the previous move was not an increase
            else
                change = change * 1.5;      %Increases the step size
                up = 1;                     %The last step was now a increase
            end

            %If the best fit is not the top row and the error changes by enough
            %the loop continues with no change in the step size
        else
            up = 0;
        end
                                            
        %When the step size is reduced to a very small amount the loop is
        %exited, for the given data the best fit has been found.
        if change < firstChange * changeCutoff
            break;
        end

        %This section outputs a current comleted percentage of the fitting,
        %it takes the percentage of lines included and then adds the
        %percentage it is between them (it know this because of the current
        %change value.
        
       
        if(up == 0)&(t0calc ==0)
            % the clc command will clear the matlab screen however will not
            % work when compiled for windows. So can be commented out when
            % compiling
            
            %clc;
             perc = ((prevLines + (log(change/firstChange)/log(changeCutoff))*(lines-prevLines))*100)/totDataSize(1);
             sprintf('%3.0f%%',perc)
        end
        
        

    end

    prevLines = lines;                  %Updates previous lines (used for percentage calculations
    lines = 2*lines;                    %Doubles the amount of lines (data) included in the fit

    %This checks if all the data was included in the previous fit, if so
    %then the best fit has been found
    if abs(lines - 2*((totDataSize(1)-1))) < 1e-5
        break;
    end

    %This checks to see if the data included in the next fit is larger than
    %80% of the complete dataset, if it is then it includes all the data,
    %the reason for this is that in the final fitting stages it takes a
    %long time each fit, so this can reduce the number by 1 in certain
    %instances.
    if  lines > totDataSize(1)*0.8
        lines = (totDataSize(1)-1);
    end


    DeltaGCat = value(1,1);
    DeltaGInact = value(1,2);
    DeltaHEq = value(1,3);
    TEq = value(1,4);
    LMS = value(1,5);

end

       
DeltaGCat = value(1,1)
if t0calc == 0
  DeltaGInact = value(1,2)
end
DeltaHEq = value(1,3)
TEq = value(1,4)

toc;

%Writes the best fit coefficients to the spreadsheet
xlswrite(filename,DeltaGCat,2,'D6');
xlswrite(filename,DeltaHEq,2,'D7');
xlswrite(filename,TEq,2,'D8');
xlswrite(filename,value(1,5),2,'D3');
if t0calc == 0
  xlswrite(filename,DeltaGInact,2,'D5');
end

%working out some statistics....NEW CODE ADDED BY KARIN and NIGEL STARTS HERE

if t0calc == 1
    model = calcT0(T,t,value(j,1),value(j,2),value(j,3),value(j,4),E0,h,Kb,R); 
else
    model = calc(T,t,value(j,1),value(j,2),value(j,3),value(j,4),E0,h,Kb,R); 
end
  totalss=sum(y.^2);                %total sum of squares
  resss=sum((y-model).^2);          %residual sum of squares
  regss=totalss-resss;              %regression sum of squares                    
  rsquare=regss./totalss;                       
  corr=mean(y.*model)./std(y)./std(model);
  
%Writes the statistics to the spreadsheet
xlswrite(filename,totalss,2,'D12');
xlswrite(filename,resss,2,'D13');
xlswrite(filename,regss,2,'D14');
xlswrite(filename,rsquare,2,'D15');
xlswrite(filename,corr,2,'D16');



%works out the confidence intervals on each of the fitted values
if highres==1
    if t0calc == 0
        button = questdlg('Calculating error statistics could take between 1 and 2 hours. Do you wish to continue?','High resolution calculation warning.','Yes','No','Yes')
        comparestring = strcmp(button,'Yes')
        if comparestring == 1
        
        
'Calculating error statistics...this could take between 1 and 2 hours'
[valueCI,error]=bootstrap(totalData,value(1,:),err,E0,Kb,R,h,webpage);

%EACAT
xlswrite(filename,error(1,1),2,'F6'); %standard deviation
xlswrite(filename,error(2,1),2,'E6'); %mean
xlswrite(filename,error(3,1),2,'G6'); %5% confidence interval
xlswrite(filename,error(4,1),2,'H6'); %95% confidence interval

%EAl

xlswrite(filename,error(1,2),2,'F5'); %standard deviation
xlswrite(filename,error(2,2),2,'E5'); %mean
xlswrite(filename,error(3,2),2,'G5'); %5% confidence interval
xlswrite(filename,error(4,2),2,'H5'); %95% confidence interval

%HEq
xlswrite(filename,error(1,3),2,'F7'); %standard deviation
xlswrite(filename,error(2,3),2,'E7'); %mean
xlswrite(filename,error(3,3),2,'G7'); %5% confidence interval
xlswrite(filename,error(4,3),2,'H7'); %95% confidence interval

%TEq
xlswrite(filename,error(1,4),2,'F8'); %standard deviation
xlswrite(filename,error(2,4),2,'E8'); %mean
xlswrite(filename,error(3,4),2,'G8'); %5% confidence interval
xlswrite(filename,error(4,4),2,'H8'); %95% confidence interval
if t0calc == 1
model = calcT0(T,t,error(2,1),error(2,2),error(2,3),error(2,4),E0,h,Kb,R); 
else
model = calc(T,t,error(2,1),error(2,2),error(2,3),error(2,4),E0,h,Kb,R); 
end

  totalss=sum(y.^2);               %total sum of squares
  resss=sum((y-model).^2);          %residual sum of squares
  regss=totalss-resss;                          %regression sum of squares                    
  rsquare=regss./totalss;                       
  corr=mean(y.*model)./std(y)./std(model);

%Writes the statistics to the spreadsheet
xlswrite(filename,totalss,2,'E12');
xlswrite(filename,resss,2,'E13');
xlswrite(filename,regss,2,'E14');
xlswrite(filename,rsquare,2,'E15');
xlswrite(filename,corr,2,'E16');
xlswrite(filename,sqrt(resss/length(y))/length(y),2,'E3'); %Standard Error
        end
    end
end

%NEW CODE STOPS HERE
        
if webpage ~= 1
    sprintf ('DeltaGCat = %20.18f',DeltaGCat)
    sprintf ('DeltaGInact = %20.18f',DeltaGInact)
    sprintf ('DeltaHEq = %20.18f',DeltaHEq)
    sprintf ('TEq = %20.18f',TEq)
    sprintf ('LMS = %20.18f',value(1,5))
end



if webpage
    
       
    %Prompt to check if the user wants to view a graph of the datapoints and
    %the best fit.
    %     graph = input('Would you like to view a graph of the fit? (Yes/No):\n','s');
    graph = 'Yes';
    Data = allData;     %Makes sure the data included in the plot is all data
    dataSize = size(Data);  %Since data has changed, size needs to be updated
    
    %Only if the user wants the graph to be displayed
    if strcmpi(graph,'Yes')
        'Plotting results'
        figure(1);
        clf;
        
      if t0calc 
            Curvefit = calcT0(Data(1,2:end),t,DeltaGCat,DeltaGInact,DeltaHEq,TEq,E0,h,Kb,R);
            plot(Data(1,2:end),Data(2,2:end));
            hold on;
            plot(Data(1,2:end),Curvefit);
            title('Fitted Curve Check','FontSize',20');
      else
          
          
        %Clears any previous data that may exist
        clear graphData;
        clear graphData2;
        count = 1;

        %This steps through all the datapoints and creates and array with the
        %x, y and z data as seperate rows.
        for i = 2:dataSize(1)
            for k = 2:dataSize(2)
                Temp = Data(1,k);
                Time = Data(i,1);
                graphData(1,count) = Temp;
                graphData(2,count) = Time;
                graphData(3,count) = Data(i,k);
                graphData(4,count) = (-calc(Temp,Time,DeltaGCat,DeltaGInact,DeltaHEq,TEq,E0,h,Kb,R)+Data(i,k)); %Finds the residual at each point
                count = count + 1;

            end
        end


        %Following line creates a meshgrid from the maximum to minimum in each
        %direction and makes 100 points in each direction
        [X,Y] = meshgrid(min(graphData(1,:)):(max(graphData(1,:))-min(graphData(1,:)))/100:max(graphData(1,:)),min(graphData(2,:)):(max(graphData(2,:))-min(graphData(2,:)))/100:max(graphData(2,:)));
        Xsize = size(X);

        %This finds the z data accoring to the best fit equation, at points
        %given y the meshgrid
        for i=1:Xsize(1)
            for j=1:Xsize(2)
                Z(i,j) = calc(X(i,j),Y(i,j),DeltaGCat,DeltaGInact,DeltaHEq,TEq,E0,h,Kb,R);
                
                %calculates the rate of reaction at each point storing it
                %in Zrate. -DH
                Zrate(i,j) = calcV(X(i,j),Y(i,j),DeltaGCat,DeltaGInact,DeltaHEq,TEq,E0,h,Kb,R);
            end
        end
        
        
     
        [XTime,YTemp] = meshgrid(Data(2:end,1),Data(1,2:end));
        XTimesize = size(XTime);
        col = 1;      
        for i=1:XTimesize(1)
            for j=1:4:XTimesize(2)
                ZAct(i,j) = calcV(YTemp(i,j),XTime(i,j),DeltaGCat,DeltaGInact,DeltaHEq,TEq,E0,h,Kb,R);
                OutPutComb(col,1) = YTemp(i,j);
                OutPutComb(col,2) = XTime(i,j);
                OutPutComb(col,3) = ZAct(i,j);
                col = col + 1;
            end
        end



        %This section outputs the fit and a residual curve to the screen.
        subplot(1,3,1);
        plot3(graphData(1,:),graphData(2,:),graphData(3,:),'.');    %Plots the datapoints
        hold on;
        mesh(X,Y,Z);                                                %Plots the equation surface in 3d
        title('Fitted Curve Check','FontSize',20');
        hidden off;
        view(145,20);

        subplot(1,3,2);
        plot3(graphData(1,:),graphData(2,:),graphData(4,:),'.');    %Plots the residual datapoints.
        title('Residual','FontSize',20);
        view(145,20);
        
        %This plots out the rates in a mesh and titles the graph-DH
        subplot(1,3,3);
        mesh(Y,X,Zrate);
        title('Rate Curve','FontSize',20);
        view(145,20);
        
        %Need to write data to sheet. So creates a sheet called Rate Data
        %and writes to that. -DH
        warning off MATLAB:xlswrite:AddSheet;

          xlswrite(filename, 'T', 'Rate Data Actual', 'A1');
          xlswrite(filename, 't', 'Rate Data Actual', 'B1');
          xlswrite(filename, 'R', 'Rate Data Actual', 'C1');
          xlswrite(filename, OutPutComb, 'Rate Data Actual', 'A2');
      end
    end
end                                 %The end of the program :D
end

